%% Contourlet变换
clc;clear;close all;

img = imread("D:\学习资料\大三下学习资料\数字图像处理\参考文献\代码部分\3.Lena图像.tif");
% figure;imshow(img)
% img = imread("D:\学习资料\大三下学习资料\数字图像处理\参考文献\代码部分\花瓶.tif");
img = im2double(img);
% NSCT 变换
y = nsctdec(img,[2,4],'pkva','9-7');   % 高斯对应2 2

% % 低频图像
% figure;imshow(y{1});
% 
% % 一尺度图像
% for i = 1:4
%     figure;imshow(y{2}{i},[]);
% end
% 
% 二尺度图像
% for i = 1:16
%     figure;imshow(y{3}{i},[]);
% end


%% 水印图像获取并进行奇异值分解
w = imread("D:\学习资料\大三下学习资料\数字图像处理\参考文献\代码部分\Cameraman.tif");
w = w(:,:,1:3);
w = rgb2gray(w);
w = im2double(w);

% w = arnold_transform(w,5);

[U_w,S_w,V_w] = svd(w);


%% 原始图像四尺度NSCT变换，在高频

t = y;

S_w1 = S_w(1:32,1:32);
S_w2 = S_w(33:64,33:64);

[U_f1,S_f1,V_f1] = svd(t{1});
[U_f2,S_f2,V_f2] = svd(t{3}{5});
%figure;imshow(y{1})

alpha = 0.1;
temp = S_f1(1:32,1:32) + 0.1 * S_w1 ;
S_f1(1:32,1:32) = temp;

temp = S_f2(1:32,1:32) + 0.2 * (S_w2) ;
S_f2(1:32,1:32) = temp;

% a = diag(S_w);
% for i = 1:8
%     for j = 1:8
%         [U_f,S_f,V_f] = svd(t{1}(i*64-63:i*64,j*64-63:j*64));
%         S_f(1,1) = S_f(1,1) + 0.2*a((i-1)*8+j);
%         t{1}(i*64-63:i*64,j*64-63:j*64) = U_f * S_f * V_f';
%     end
% end


%% 图像复原 逆NSCT
t{1} = U_f1 * S_f1 * V_f1';
t{3}{5} = U_f2 * S_f2 * V_f2';
f_watered = nsctrec(t,'pkva','9-7'); % f_watered即为添加水印后的图像

% 图像重建
figure;imshow(f_watered,[]);
% figure;imshow(img - f_watered,[]);

% % 测试
% y_watered = nsctdec(f_watered,[2,3],'pkva','9-7');
% [U_watered,S_watered,V_watered] = svd(y_watered{3}{5});



%% 水印图像不可见性测试，PSNR
f1 = double(img);
f2 = double(f_watered);

% 计算差值
diff_f = f1 - f2;

% 计算差值的平方
diffImageSquared = diff_f .^ 2;

% 计算MSE：差值平方和的均值
mse_value = sum(diffImageSquared(:)) / numel(f1);

psnr_value = 0;
if mse_value == 0
    psnr_value = Inf;
else
    psnr_value = 10 * log10((1^2) / mse_value);
end

disp(['PSNR: ', num2str(psnr_value), ' dB']);
if psnr_value > 40
    disp("水印图像不可感知性良好");
else
    disp("水印图像不可感知性较差");
end


%% 无攻击提取水印测试
NC_JPEG_compressed = water_extracted(img,f_watered,w)


%% 最佳嵌入强度选取
% [PSNR,NC] = myfun_best_alpha(img,w);
% myfun_plot_alpha(PSNR,NC);

%% 鲁棒性测试
% JPEG压缩测试
Q = 80; % 质量因子=80
imwrite((f_watered),"Q=80.jpg",'Quality',Q);
f_JPEG_compressed = imread("Q=80.jpg");
figure;imshow(f_JPEG_compressed,[]);
%f_JPEG_compressed = rgb2gray(f_JPEG_compressed);
NC_JPEG_compressed = water_extracted(img,(f_JPEG_compressed),w)

% JPEG压缩测试
Q = 50; % 质量因子=80
imwrite((f_watered),"Q=50.jpg",'Quality',Q);
f_JPEG_compressed = imread("Q=50.jpg");
figure;imshow(f_JPEG_compressed,[]);
%f_JPEG_compressed = rgb2gray(f_JPEG_compressed);
NC_JPEG_compressed = water_extracted(img,(f_JPEG_compressed),w)

Q = 20; % 质量因子=20
imwrite((f_watered),"Q=20.jpg",'Quality',Q);
f_JPEG_compressed  = imread("Q=20.jpg");
figure;imshow(f_JPEG_compressed,[]);
NC_JPEG_compressed = water_extracted(img,f_JPEG_compressed,w)

% 中值滤波
f_median_fliter = medfilt2(f_watered,[2,2]);
figure;imshow(f_median_fliter,[]);
NC_median_fliter = water_extracted(img,f_median_fliter,w)

% 中值滤波
f_median_fliter = medfilt2(f_watered,[3,3]);
figure;imshow(f_median_fliter,[]);
NC_median_fliter = water_extracted(img,f_median_fliter,w)

% 均值滤波
f_amf = imfilter(f_watered,fspecial("average",[2,2]));
figure;imshow(f_amf,[]);
NC_amf = water_extracted(img,f_amf,w)

% 均值滤波
f_amf = imfilter(f_watered,fspecial("average",[3,3]));
figure;imshow(f_amf,[]);
NC_amf = water_extracted(img,f_amf,w)

% 高斯噪声
noise = imnoise2('gaussian',512,512,0,0.005);
f_gauss_noise = (f_watered) + noise;
figure;imshow(f_gauss_noise,[]);
NC_gauss_noise = water_extracted(img,f_gauss_noise,w)

% 高斯噪声
noise = imnoise2('gaussian',512,512,0,0.01);
f_gauss_noise = (f_watered) + noise;
figure;imshow(f_gauss_noise,[]);
NC_gauss_noise = water_extracted(img,f_gauss_noise,w)

% 高斯噪声
noise = imnoise2('gaussian',512,512,0,0.02);
f_gauss_noise = (f_watered) + noise;
figure;imshow(f_gauss_noise,[]);
NC_gauss_noise = water_extracted(img,f_gauss_noise,w)

% 指数噪声
noise = imnoise2('exponential',512,512,25,2);
f_exp = double(f_watered) + noise;
figure;imshow(f_exp,[]);
NC_s_p = water_extracted(img,f_exp,w)

% 椒盐噪声
noise = imnoise2('salt & pepper',512,512,0.02,0.02);
f_s_p = double(f_watered) + noise;
figure;imshow(f_s_p,[]);
NC_s_p = water_extracted(img,f_s_p,w)

% 直方图均衡攻击
f_histogram_equal = imadjust(f_watered);
figure;imshow(f_histogram_equal,[]);
NC_histogram_equal = water_extracted(img,f_histogram_equal,w)

% 旋转攻击
f_rotated = imrotate(f_watered,5,'bicubic', 'crop');
figure;imshow(f_rotated,[]);
NC_rotated = water_extracted(img,f_rotated,w)

% 裁剪攻击
f_cutted = f_watered;
f_cutted(1:128,1:128) = 0.5;
figure;imshow(f_cutted,[]);
NC_cutted = water_extracted(img,f_cutted,w)

% % 颜色污染
% f_stain = f_watered + 0.1;
% figure;imshow(f_stain,[]);
% NC_cutted = water_extracted(img,f_stain,w)

% 平移操作
B = imtranslate(f_watered,[256,20]);
B(1:20,:)=f_watered(512-19:512,:);
B(:,1:10)=f_watered(:,512-9:512);
figure;imshow(B,[]);
NC_cutted = water_extracted(img,B,w)

% 
f_tran = imtranslate(f_watered,[50,50]);
f_tran = f_tran(51:512,51:512);

% 缩放攻击
scaleFactor = 0.9;
f_scale = imresize(f_watered,scaleFactor);
figure;imshow(f_scale,[]);
NC_cutted = water_extracted(img,f_scale,w)


